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Abstract:  Formulation  is  given  for  efficient  parabolic  equation  solution  of  radiowave 
propagation  in  inhomogeneous  atmosphere  and  over  irregular  terrain.  Both  standard 
and  wide  angle  parabolic  equation  derivations  are  presented.  Impedance  boundary 
conditions  are  used  to  characterize  the  ground.  A  tropospheric  boundary  condition 
based  on  the  exact  solution  of  Schrodinger  equation  in  a  quarter  plane  is  derived.  To 
permit  efficient  modeling  of  the  irregular  boundary,  the  parabolic  equation  together 
with  the  boundary  conditions  are  transformed  into  a  numerically  generated  curvilinear 
coordinate  system.  Finally,  formulation  is  presented  for  a  finite  difference  solution 
using  Crank-Nicolson  implicit  scheme. 


1.      Introduction 

It  is  well  known  that  multipath  fading  can  significantly  effect  the  link  reliability  in  a 
communications  system  or  target  detectability  in  the  case  of  radar.  The  path  between 
a  transmitted  and  a  receiver  is  often  obstructed  by  natrual  or  man-made  obstacles 
such  as  hills,  buildings,  atmospheric  layers,  trees,  rain,  fog,  etc.  Propagation  outage 
due  to  multipath  fading  depends  in  a  complicated  manner  on  propagation  climate, 
terrain  features,  path  length,  radio  frequency,  and  fade  margin.  In  the  case  of  atmo- 
spheric multipath  fading,  interference  due  to  two  or  more  super-refracted  rays  arriving 
at  the  receiver  via  different  paths  can  lead  to  a  complete  loss  of  signal.  Other  patho- 
logical phenomenon  such  as  obstruction  fading  (caused  by  sub-refractive  atmospheric 
effects)  and  ducting  (caused  by  extreme  super-refractive  effects)  are  also  possible. 
Such  phenomenon  are  more  common  in  warm  and  tropical  climates,  particularly  near 
shores,  where  elevated  inversions  are  formed  easily  due  to  the  large  temperature  and 
partial  pressure  differentials.  Reflection  multipath  fading,  which  is  due  to  interfer- 
ence between  the  direct  and  the  ground  reflected  ray  depends  strongly  on  the  terrain 
geometry  and  ground  constants.  Moreover,  elevated  terrain  features  could  completely 
mask  a  receiver  from  a  transmitter  leading  to  severe  loss  of  signal  (in  some  cases,  it 
is  advantageous  to  site  antennas  behind  hills  to  provide  shielding  against  undesirable 
interference).  It  is  very  important  to  assess  the  effects  of  environment  on  the  link.  A 
computer  model  that  can  take  into  account  a  given  refractive  index  profile,  terrain 
elevation  data,  and  varying  ground  parameters  will  be  very  helpful  in  predicting  the 
link  performance. 

In  this  report  we  present  formulation  details  for  an  efficient  numerical  solution  of 
wave  propagation  in  an  inhomogeneous  atmosphere  and  over  irregular  terrain  using 
parabolic  equation. 

Unlike  all  other  previous  formulations  of  the  parabolic  equation,  we  will  use  a  modified 
Helmholz  equation  for  propagation  in  an  inhomogeneous  atmosphere  as  suggested  by 
Maxwell's  equations  (all  previous  formulations  use  a  Helmholtz  equation  which  is 
only  true  for  fields  in  a  homogeneous  medium).  Because  the  parabolic  equation  is  a 
full-wave  method,  it  will  include  all  aspects  of  wave  propagation  such  as  reflection, 
refraction,  diffraction,  and  surface  wave  propagation.  In  this  respect  it  is  far  superior 
to  the  commonly  used  ray  method. 

Parabolic  equation  approximation  to  an  elliptic  partial  differential  equation,  which 
the  true  fields  satisfy,  has  proven  to  be  a  viable  approach  for  studying  propagation 
problems  in  underwater  acoustics.  The  method  is  just  gaining  popularity  with  the 
electromagnetic  community.  Although  the  parabolic  equation  regards  waves  as  es- 
sentially traveling  one-way,  it  allows  a  rapid  solution  of  the  fields  by  way  of  marching 
along  the  range  starting  from  an  initial  range.  Another  advantage  of  the  PE  method 


compared  to  the  ray  methods  is  that  it  is  valid  even  in  the  shadow  region  where  the 
simple  ray  methods  completely  break  down.  Furthermore,  it  appears  to  be  the  only 
practical  method  for  predicting  propagation  over  long  ranges  (greater  than  1  km)  over 
a  wideband  from  HF  (a  few  MHz)  through  SHF  (a  few  tens  of  GHz).  The  method  is, 
however,  not  without  limitations.  In  its  standard  form,  the  accuracy  of  the  method 
is  limited  to  waves  traveling  essentially  within  ±10°  from  horizontal.  Furthermore, 
treatment  of  the  boundary  conditions  on  the  uneven  terrain  is  difficult. 

The  method  we  propose  to  use  will  attempt  to  remove  both  of  these  defficiencies. 
Firstly,  we  use  a  Helmholtz-like  elliptic  equation  to  describe  the  fields  and  arrive 
at  a  wide-angle  parabolic  equation  subject  to  certain  approximations.  To  facilitate 
imposition  of  the  boundary  conditions  on  the  irregular  terrain,  the  equations  will  be 
transformed  to  a  body-fitted  curvilinear  coordinate  system.  The  PE  will  be  solved 
using  finite  differences  on  a  non-rectangular  mesh.  This  is  a  major  departure  from 
previous  approaches  which  will  not  only  make  the  method  more  efficient  but  also 
more  accurate. 

In  Section  2,  we  derive  the  exact  equations  satisfied  by  2D  fields  in  an  inhomogeneous 
atmosphere.  Impedance  boundary  conditions  are  used  to  characterize  the  ground.  In 
Section  3,  we  present  details  on  the  impedance  boundary  conditions.  Starting  from 
the  exact  equations  presented  in  Section  2  for  the  fields,  we  derive,  in  Section  4,  a 
parabolic  equation  (PE)  valid  for  narrow  angle  propagation.  This  case  is  termed  as 
the  standard  PE.  The  standard  PE  is  generally  valid  for  propagation  angles  that 
are  within  ±10°  from  horizontal.  To  accomodate  waves  traveling  at  larger  angles,  we 
present  the  derivation  of  a  wide  angle  PE  in  Section  5.  To  truncate  the  computational 
domain  we  derive  boundary  conditions  on  an  upper  boundary,  which  are  termed  as 
the  tropospheric  boundary  conditions.  The  derivation  is  based  on  the  solution  of 
Schrodinger  type  parabolic  equation  in  a  quarter  plane  x  >  0.  y  >  0.  This  is  presented 
in  Section  6. 

For  an  efficient  numerical  implementation,  we  transform  the  differential  equation  and 
boundary  conditions  to  a  curvilinear  coordinate  system.  This  is  presented  in  Section 
7.  Details  on  the  numerical  generation  of  the  curvilinear  coordinate  system  are  given 
in  Section  8.  Finally,  in  Section  9.  we  present  steps  leading  to  a  finite  difference 
solution  of  the  equations  using  a  Crank- Nicholson  implicit  scheme. 


2.      Solution  of  2D  Fields  in  an  Inhomogeneous  Medium 

Consider  an  electric  source  producing  fields  in  an  inhomogeneous  region  as  shown 
in  the  figure  below.  Let  us  assume  that  both  the  sources  and  the  medium  are  two- 
dimensional  in  nature  in  that  all  quantities  are  independent  of  the  z-coordinate.  As 
in  the  case  of  a  homogeneous  medium  the  fields  can  be  decomposed  into  a  TEZ  case 
(vertical  polarization)  and  a  TMZ  case  (horizontal  polarization).  It  is  assumed  that 
propagation  takes  place  in  the  rry-plane. 
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From  Maxwell's  equations,  we  have 

V  x  E    =    —ju[iH 

V  x  H    =    jueE  +  J 


(1) 
(2) 


In  the  case  of  vertical  polarization^  the  fields  could  be  written  in  terms  of  the  z- 
component  of  the  magnetic  field,  11  —  zll2\  (TEZ  fields).  Substituting  into  (2)  we 
have 


V  x  H  =  V  x  (£//,)  =  Vllz  x  z  =  jucE  +  J  ==>  eE  = 
In  a  source-free  environment,  we  have 


V//z  xz-J 


tE  = 


VHZ  x  z 


Substituting  into  (1)  we  get 


V  x  E  =  —V  x x  z)  =  -—V  • =  -wull 


The  equation  satisfied  by  the  magnetic  field  is  then 

y  — —  j  =  -jlvhH2 


(^]+u-V^  = 


Vertical  Polarization 


(3) 


Once  the  magnetic  field  is  determined  the  electric  field  is  given  by 


E 


jut 


(4) 


Note  that  IL  does  not  satisfy  the  Helmholtz  equation  unless  c  is  constant. 
For  horizontal  polarization  on  a  similar  analysis  with  E  —  zEz  shows  that 

Horizontal  Polarization 


(5) 


If  the  medium  is  non-magnetic,  fi  —  /v0  and  Ez  satisfies  the  Helmholtz  equation.  We 
will  characterize  the  ground  in  terms  of  impedance  boundary  conditions  [3]. 

We  may  combine  the  vertical  and   horizontal  cases  shown  in   (3)   and   (5)  into  an 
equation  of  the  form 

V  •  (aW)  +  3v  =  0  (6) 


where 


1 
e 


a     —     < 


TE  or  Vertical  Pol. 
TM  or  Horizontal  Pol. 


0    -     QA-oV(.r.y)>0 
H2        TE  Pol. 

TM  Pol. 


(7) 

(8) 
(9) 


The  quantity  n(x,y)  is  a  position  dependent  refractive  index  of  the  medium.    The 
partial  differential  equation  (PDE)  given  in  (6)  is  elliptic,  for  we  have 

dx  \    dx)      dy  \    dy) 


or 


Q 


0V      d\> 
+ 


dv  da       dij'  da 
ox   ox       dy    dy 


dx2       dy2 

The  discriminant  for  the  above  PDE  is  —1  <  0  implying  elliptic  nature  of  the  equation. 
Equation  (6)  could  be  expressed  as 


1    da  1    da 


a   ox  a   oy 


a 


For  a  non-magnetic,  loss-less  medium,  we  have 
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tot, 


Vertical  Pol. 


a  =  < 


1 


—  ,         Horizontal  Pol. 
^     Co 


so  that  for  vertical  polarization 


a   dx 


d    (\ 


Zt 


dx  \tr 

2   du_ 
n   dx 


1    dcT  1     dn2 

cr    dx  n2    dx 

(n  is  the  refractive  index' 


Similarly, 
Letting 


1    da  2   dn 

a   dy  n   dy 


a^x.y) 


^(x-.y)   = 


2  dn 

77     dx 

0 
2     #77 

n   dy 
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Vertical  Pol. 
Horizontal  Pol. 
Vertical  Pol. 

Horizontal  Pol. 


equation  (10)  may  be  expressed  in  the  form 

V2i/>  +  a\4>x  +  a2ipy  +  fcJnV  =  0 


(10) 


(11) 


3.      Impedance  Boundary  Condition 

Impedance  boundary  condition  relates  the  tangential  components  of  electric  and 
magnetic  fields  at  the  interface  of  two  media.  If  v  is  a  unit  normal  and  s  is  a  unit 
tangent  as  shown  in  the  figure  below,  the  boundary  conditions  are  given  by  [3] 


s 


r 


s 


\ 


V 


s  —  i  cos  9  -f  y  sin  9 
v  =  —x  sin  9  +  y  cos  9 

x  =  s  cos  9  —  0  sin  9 
y  =  s  cos  9  -\-  v  cos  9 


0  x  (0  x  £)  =  -t]0Asu  x  H  (12) 

where  As  =  Zs/Vo  's  ^nc  surface  impedance  normalized  to  the  free  space  value  r)0. 
The  equation  may  also  be  expressed  as 


(0  ■  E)v  —  E  =  -T]0AS0  x  II  =>  0  x  E    =    tj0AsO  x  (0  x  H) 

=    VoX  \{0-H)0-H 


(13) 


The  surface  impedance  is  determined  from  the  intrinsic  impedance  of  the  medium 
by  considering  plane  wave  reflections  from  the  interface.  The  complex  propagation 
constants.  7^  72,  and  the  intrinsic  impedances  771  and  772  in  terms  of  the  media 
constants  are  indicated  in  the  figure. 
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eiiCo^vi^,  rt     72  =  -k&r2erc2 


— 7~ =  Vo\  , 


V2     =     Vo< 


f/fr2 
Cr2 
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According  to  Snell's  law,  we  have  7j  sin  9l  =  *y2  sm  &t 


The  plane  wave  reflection  coefficients  for  the  vertical  and  horizontal  polarizations,  Ry 
and  Rh  are  [7] 
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For  the  special  case  of  )i\  =  //2  =  ^0-  Ci  =  0, 
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(15) 


Crl 


a:  =  ,/ 

V  Cr2  -  JCTr2 

For  normal  incidence  v,  =  90° 

a?  =  a;'  - 

V  e^2  -  J0>2 

For  the  2-D  case,  impedance  boundary  conditions  for  vertical  polarization  can  be 
simplified  as 


Taking  a  dot  product  on  both  sides  with  z 
z-{uxE)  =  -Vo^]'Hz 

Since  z  x  i>  =  —5,  we  have 
-s-E=  -riots*  Hx 

Substituting  from  equation  (4)  for  E,  we  get 

.    z  x  VHZ 

s  •  : = 

JUit 


■voa:hz 


u-VH2  =  jcoerjoAl'H,  =  jk0trA}sHz 
i.e.,  the  above  can  be  put  conveniently  in  the  form 


dH2 
dv 


]h(rAvsHz  =  0 


IBC  Vertical  Pol. 


A  similar  analysis  for  horizontal  polarization  yields 


dE>        i         l    r       n 


IBC  Horizontal  Pol. 


This  may  also  be  obtained  by  resorting  to  duality. 
For  a  perfectly  conducting  material,  we  have  A}J'^    =  0  and 

dll. 


=    0 
dv 

E,    =    0 


Vertical  Pol. 
Horizontal  Pol. 


(16) 


(17) 


4.      Standard  PE  Derivation 

We  will  make  some  approximations  and  cast  (11)  in  the  form  of  a  parabolic  equa- 
tion which  permits  a  rapid  numerical  solution. 

e~jk0x 

Let  xl){x,y)  =  — -=-u{x,y) 


Then 
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%     =     w 


e-J'*°x 


I 

•        Wyx   =  ^rj/   ~   {Uly  ~  jkQUy)- 


X   — ->   OC 


3—jkox 

y/x 
P-jk0x 


Vyy        —        U, 


Substituting  into  (11)  we  get 


•      uXT  ~  (uIT -2jk0uT  -  k0u)- 


^ 


Ku  +  Uyh  -  2jk0uT  +  G!i/r  -  2jk0axu  +  a2uy  +  (n2  -  l)fcju  =  0 


uxx  +  ^yV  +  («i  -  2jk0)ux  +  a2i/y  -f  in2  -  1  -  2j  —  j  k\u  -  0 
If  now  we  impose  the  approximation  that 

/     9  o  \  1  /  2 


we  obtain 


"  (ai-2j*b) 


2  ~ ]  _2jV")  ^°w 


or 


(18) 


This  is  the  exact  form  of  narrow  angle  PE  approximation.    We  would  also  like  to 
express  the  impedance  boundary  condition  in  terms  of  the  'u'  functions. 


x    =  s  cos  9  —  v  sin  6 
dx 

Xv  =  7—      =  v  ■  Vx  =  V  ■  X 
OV 

=  —  sin# 

v    =  —  x  sin  6  +  y  cos  6 


For  vertical  polarization  we  have 
OIL 


dv 


-jkoeTAlH;  =  0 


dH2 
dv 


jk0n2A]'Hz  =  0 


Hz(x.y)     = 


,-jkox 
y/x 


■u(x.y) 


dllz  (  uxu   \  e 

=     \uu  -  jk0xuu  - 


-jkox 


dv 


(uu  -  jk0xuu] 


2x3/2, 

,-jk0x 


x  — ►  oc 


u»  -  jk0{n2As  +  xv)u  =  0 


\/x 

IBC  Vert.  Pol. 


Similarl\'  for  horizontal  polarization,  we  have 


IBC  Horz.  Pol. 


We  combine  the  two  by  defining 


Cl 


■jk0{n2A]'  -  sinfl)        Vert. 


-jfco 


A" 


sin  0  Horz. 


and  writing  as 


uL,  +  c-iu  =  0 


(19; 


The  parabolic  equation  given  in  (18)  is  valid  for  propagating  angles  close  to  horizon- 
tal (±10°  in  practice)  [1].  To  accomodate  waves  at  higher  angles  we  would  need  a 
wide  angle  parabolic  equation  whose  derivation  is  accomplished  through  a  pseudo- 
differential  operator  formalism  [1]. 
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5.     Wide  Angle  PE  Derivation 


Let  us  assume  that  the  media  constants  are  independent  of  horizontal  range  so  that 
a(x,y)  =  a(y).  n{x,y)  =  n(y).  We  then  have  from  (17) 

uxx  +  uyy  -  j2k0ux  +  -  —uy  +  (n2(y)  -  tfkZu  -  0 

a   dy 


Let 


OX 


\ 


1  d2 


+ 


1    da  d 


+  n- 


Pseudo  differential 
operator  in  y 


(20) 


kl  dy2       kla  dy  dy 

Using  this  notation,  the  above  PDE  ma}-  be  expressed  as 

P2  -  2jk0P  +  (Q2  -  l)kl]  u  =  0 

which  we  may  factorize  as 

(P  -  jko  -  jk0Q)(P  -  jk0  +  jkQQ)u  =  0 

To  see  this  we  expand  the  operator  on  the  left  hand  side  to  get 

P2  -  jk0P  -  jkoQP  -  jk0P  -  k20  -  k20Q 
+  jk0PQ  +  k20Q2  +  k2QQ 

Since  PQ  —  QP,  we  have  the  desired  result.  The  first  operator  in  (20)  denotes 
an  incoming  wave  (w.r.t.  x)  and  the  second  an  outgoing  wave.  We  retain  only  the 
outgoing  wave  to  obtain 

Pu  =  -jk0{Q-l)u  (21) 

The  square  root  operator  Q  is  global  in  nature  and  we  would  like  to  make  some 
approximations  to  derive  a  local  operator  from  it.  (It  is  global  because  when  expanded 
in  terms  of  series,  it  will  contain  terms  of  all  derivatives).  Now 


Q 


k\  a   dy   dy       k%  dy2  _ 


1/2 


Let  us  rewrite  Q  as 


Q  = 


l  +  ["2(y)-l]  +  - 


1       da         d 


+ 


d2 


<C1  normaLlv 


a  d(k0y)d(k0y)       d(k0y) 


iall 


1/2 
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which  may  be  expressed  as 
where 


g  =  vT+v 


V    =    n2  -  1  + 


1     da    d        1     cP 

+ 


k2a   dy   dy       k20   dy2 
=     a  small  operator! 

Treating  V  as  an  algebraic  factor  <  1.  we  may  derive  the  following  rational  approxi- 
mation (pade(] .  1 )) 

Q  =  >/TTT    %     |—      (Claerbout) 

4  +  3V  2V 

=  JTV  =  1  + 

2V 


4  +  V 


so  that 


Q-\  = 


4  +  V 


Substituting  into  (21)  we  arrive  at 


Pu  =jkQ{Q-  l)u=»  Pu  =  - 


4  +  V 


or 


4  +  V)Pu  =  -2jk0Vu 


1    da 


Using  the  fact  that  a2  — — .  we  get 

q    dy 


0         d2 

V  +  3^  +  02c^+cV 


ux  =  -2jfc0 


5         0s 


V-l)^  +  a2—  + 


dy      dy2 


(22) 


The  above  equation  is  a  wide-angle  parabolic  equation  valid  for  propagation  up  to 
±20°  [1].  Other  approximations  could  be  obtained  by  considering  higher  order  padt 
approximants. 
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6.      Boundary  Condition  on  the  Upper  Boundary 


To  truncate  the  computational  domain,  we  consider  a  point  high  enough  where  the 
atmosphere  is  homogeneous  with  n  —  \.  The  governing  equation  in  the  homogeneous 
region  becomes 

3 


U-r  =   - 


2fcn    yV 


Let  us  derive  boundary  conditions  on  a  horizontal  interface  y  =  yo.  For  the  sake 
of  simplicity  we  will  derive  boundary  conditions  of  the  mixed  type  on  the  interface 
y  —  0  (instead  of  y  —  yo)  given  the  initial  data  on  the  line  x  =  0  and  boundary  data 
on  the  line  y  =  0.  Our  derivation  is  based  on  the  use  of  Fourier  sine  transforms  as 
suggested  in  [2].  Although  the  basic  philosophy  of  our  approach  coincides  with  that 
in  [4],  some  of  the  details  and  the  final  results  are  slightly  different  from  the  latter. 


Consider  the  parabolic  equation  uyy 
y  >  0,  where  k  is  a  complex  constant, 


2jkux  =  0  in  a  homogeneous  region  x  >  0, 


y 


u  =  fiv) 


U 


vy 


=    zjk  U. 


/ 


-*- 


subject  to  the  initial  condition  u(0,  y)  =  f{y).  0  <  y  <  oo,  and  the  boundary  condition 
u(x.Q)  =  g(x).  0  <  x  <  oc.  We  assume  that  i/(x,oo)  — >  0  and  uy(x,oc)  — ►  0.  The 
equation 

uyy  =  2jkux  (23) 

is  of  Schrodinger's  type.  We  will  treat  the  lossless  case  having  a  real  value  of  k  as 
the  limiting  case  of  the  lossy  problem  having  k  =  k0  —  je,  e  >  0.  We  will  solve  the 
problem  using  Fourier  sine  integral. 


Let 


Us(x,\)  =  J-         u(x,y)s'm\ydy 


(24) 


14 


Then  using  integration  by  parts,  we  see  that 

(2    f°°  [2   f 

-/      uyy(x,y)smXydy    =     \  -  <uy(x,y)  sin  \y 

77  Jo  V   7T 


V=0 


roc 

-A2y      u(x.y)smXydy 


—  Xu(x.  y)  cos  Xy 
Because  uy(x,oc)  — >  0  and  u(x,oc)  — ►  0.  we  have 

\j- J     yyy{x.y)su\  Xudy     =     J -Xu(x.O)  -  X2l\{x,  X) 


■Xg(x)-XH\(x.X) 


(25) 


Multiplying  both  sides  of  (23)  with  Wz/tt  sin  Ay.  integrating  over  y  =  0  to  oc  and 
making  use  of  (25),  we  have 


or 


lhg(x)-X2Us(x.X)  =  2jk^-i-s(x,X) 


9    XT  I  V.  ^    XT,         XN  A         /2      ,      v 


We  may  rewrite  the  above  equation  as 
5 


dx 


C/s(2:,A)e(A2/2^ 


^i-w^^ 


(26) 


Replacing  the  dummy  variable  x  in  (26)  with  r  and  integrating  both  sides  over  7  =  0 
to  x,  we  arrive  at 


Us{T,X)e^l2^ 


X       2 


or 


But 


17.(*,A)  =  C/s(0,A)e-^/2^  +  i-i/*   T    <,(r)e^^— W 

2jk\-k  Jt=o 

£4(0,  A)     =     J-J     u{0,y)smXydy 

=     \~  r  f(y)smXydy  =  Fs(X) 

\    IT  JO 
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U,(x,  A)  =  F,(\)e-<*'W  +  ±Jlr    9(r)eWW-*dT  (27) 

Zj  K  V    7T  ./t=0 


Finally,  taking  the  inverse  sine  transform  on  (27)  we  get 
u(x,y)     =     J-         sm\yUs(x,\)d\ 

=     ^  T  Fs(X)e-^/2^sm\yd\ 

V    7T  JX=0 

+-4t  r  A  sin  Ay  f   ^(T)c<Aa^'*><T-*)rfTdA 

7r  2jA:  Ja=o  ./t=o 

=     -/        /      /(r)sinAr^/re-(A2/2^)lsinAydA 

/*x  /*oo 

/      y(r)/"     AsinAj/^'/^-^Adr 

7T  ZJK  Jt  =  0  J\  =  o 


2    1 

+- 


=     -  r    f(T)  r    [cos{t  -  y)\  -  cos{t  +  y)\]e-{x2/2jk)Td\dT 
■k  Jt=0  Jx-o 

+  1  f    9{r)\  (-4)  f^  cos\ye^2W-*)d\dT 

7T  Jt  =  0  jk    \      oy  J   JX-Q 

roc  roc  2 

=     -  /       f(T)  /       [cos(?/  -  r)A  -  cos(y  +  r)A]  e"(A  /2j*)r<A  dr 

-A:  f  ^(-)|-{rcos^e"(A2/2j/c)(x"T)^}^ 

Defining 

TOO 

A'(x.y;  x0.y0)  =   /      cos{y  -  y0)\e-{x  /2jk){r-To)d\.      for  x  >  x0. 
Jo 

we  write  the  expression  for  u(x,y)  as 

1    f°° 
u{x,y)     =     -  /(r)[A'(x,y;  0,r)  -  A'(x,y;  0,-r)]</r 

— yf  ^Wfl-^(*.y;T.o)rfT  (28) 

7cjk  Jt=o         oy 
We  now  evaluate  the  integral  for  K.  Consider 

roo  . 

1(a)    =      /      cosaAe-(A /2jJ:>(x-I0^A,  i  >  x0 

Jx=o 
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=      f°°  cosaXe^^^'^UX 
J\=o 

J\-0 


Because  of  the  exponential  decay,  we  may  differentiate  under  the  integral  sign  to 
obtain 


d  r^ 

da  J\=<. 


XswaXe-^^-X0^k^-^d\ 


—      I      sin 
Jx=o 


qA 


d_ 


:-(A2(x-i0)/2|fc|2)(f-jfco) 


(x-x0)(c-jk0)/\k\* 


]-dX 


(x  -  x0){i  -  jk0)/\k\2 


a\k\ 


,=0      (x-x0){e-jko) 


rcosa\e-^x-x°*-jk°M2Wd\ 

Jo 


ja\k 


(x 


\i2L-r 

x0)km  Jo 


cosaye-x*iT-xo)/2>kd\ 


jak 


(i  -  x0) 


I(a 


dl(a)          jak      ,.    . 
-P  +  -^ r/(a)  =  0 

da  [x  —  Xq) 


or 


la 


/(a)e^^/(2(x-xo))j   =  o  =»  /(tt)  =  /(tt  =  0)e-Jo^/(2(x-xo)) 


Now  7(o  =  0)  can  be  written  as 

Jo 
Let  us  view  this  integral  in  the  complex  A-plane. 


Complex  A-Plane 


(29) 
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For  the  integral  to  converge  for  (x  —  x0)  >  0,  Re[X2(e  —  jk0)}  >  0,  i.e., 

Re[\2{-jk*)]  >  0  =>  Irn{\2k*)  >  0 


or 


or 


or 


Now 


0  <  Arg(A2r)  <  it 


0  <2Arg(A)-  Arg(£)  <  ?r 


-Arg(fc)<Arg(A)<  |  +  -Arg(A:) 


Arg(fr)  =  -tan"1  (  —  )  %  -—  .  -«1 


2Av 


<  Arg(A)  <  £  - 


2       2Jt0 
7(q  =  0)  =  /  e"A2(— o)(-i**)/2|*PrfA 

where  C  lies  in  the  region  of  convergence. 

In  particular,  let  us  choose  a  line  from  0  to  oc  along  the  line  Cw  defined  by 


A  = 


1   2jk 

X  —  Xq 


-v         w  =  0  to  oc 


X-Wv(X)  A 


TV  € 


(30) 


(31 


Re  (A) 


On  this  path, 


7(a  =  0) 


Arg(A)  =  I  +  -ArgW  =  --- 


2jk 


(* 


—  Xq)   Jw=0 


e  w  dw 


2(x  -  x0) 
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1(a) 


xjk 


2(x  -  xq) 


■jo7k/(2(x-x0)) 


K(x,y;  x0,yQ)  =    U^e-^^m^o)) 

V  2{x  -  x0) 


(32) 


Substituting  this  into  (28).  the  field  u(x,y)  for  the  initial-boundary  value  problem  is 
given  by 


'vjkJT^og^dy\} 


-JHy+T)2/{2x)\dT 

e-jky*/(2(x-T))\   dr 


t  '  '  "  —  '  —  C 

■Kjk 


It  is  easy  to  see  that 
d2K 


dy< 


and,  so  for  x  ^  Xq 


r  -A2  cos[(y  -  y0)X}t-x2{l-T^2^d\ 
Jo 


d2K  dK  dK_        1    d2K 

dy2  dx  dx        2jk  dy2 


which  is  the  same  equation  satisfied  by  v.  Now 

c-jk(y-r)7/(2x)   _   e-jk(y  +  r)?/(2x) 


<*>v)  =  \i^z  r./w 


2ttx  Jt=o 
1 

7TJ 


d 


l      fx  d 
y  /      g(r)  —  K(x.y:   r.Odr 

7T7V  Jt  =  0  Oil 


du 
dy 


{x-.y) 


jk    r  n-2jkT, 


y^0+ 


2TX  Jo 

1 


f(T)t^Le-Jk^/(2X)dT 


d2 


I      f1  o- 

T  /       9{T)jr^K{x.y,   r,0) 

Trjk  Jt-o  oyz 


dr 


y— 0+ 


M-  r  f[T)2le-^l^dr 
V  2kx  Jo     JK    '   dr 

■—  f g{r)2jk—K{x,y]  r,0) 

7TJK  JO  OX 


dr 


y-+  0+ 


(33) 


(34) 
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in  view  of  (34).  Furthermore,  we  note  that 


—  K(x,y;x0,yo)  =  -^—K(x0y\x0y0] 

OX  OXq 


Using  this  and  evaluating  the  first  integral  by  parts,  we  get 

-jH/T(r)e-^2/(2*)rfT 


oy 


2J^{fir)e-^2l^ 


v-o+ 


7TX 


r=0 


d 


i  rx  a 

+-  /      g{T)—K{x,y-  r,0 

7T  Jt=Q  OT 


dr 


y— 0  + 


<2jk 


77X     I  JO  J 


+  ■ 


g{T)K(x,y\T,0) 


t_0   -  f    gT(r)K(x.O:  r.0)d' 

i  JT=0 

/  — 0  + 


V     7TX  V     XX    JO 

--g(oy 


Irjk  2jk    r*      gT{r) 


I2jk 


TTX 


2x 


UW-g(O)} 


■K      Jt=0  \/x  —  T 


dr 


+  1 


WLirMe-MMdT 


T/o 


y/x 


r^4 

70      y/X  —  T         J 


From  the  compatibility  conditions  on  the  initial  and  boundary  values,  we  have 


limu(i,O)  =  0(O)  =  limw(O,y)  =  /(O) 

x— *0  y— *0 


Therefore 


du 
dy-{X'y) 


'2jk 


y-.0+ 


1*     J-  r  fr{r)e-^'l^dT  -  f    ^Lir 

7T        [\/X  Jt-0  Jt=0  \/X  —  T 


(35) 


It  is  easy  to  note  from  (35)  that  for  k  =  k0  —  jt, 


v     du(       \ 
lim  — (x,y) 


=  0. 


v^o+ 
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No  matter  how  small  the  t  is,  we  will  use  the  same  equality  for  the  limiting  case 
of  t  — ►  0.  We  will  evaluate  the  integrals  above  approximately  by  replacing  the 
derivatives  with  differences. 


Ay 


I.u; 


.,  U 


u 


u 


u       u  r  * 


u 


V»  *K  */> 


We  now  consider  the  evaluation  of  du/dy\y=o  at  x  =  xp_-[/2-  Consider  initial  data 
on  the  line  x  =  0.  Let  us  assume  that  this  initial  data  is  known  on  a  uniform  grid 
yVl  =  mAi/.  m  =  0.  1.  •  •  •.  We  approximate  the  derivative  in  the  interval  {ym-i.ym) 
bv  the  forward  difference  formula 


df(y] 

dy 


Um-\ 


±y 


y   G    (j/m-l,J/m) 


where  um  =  u(0,ym). 
Then 


Let 


1 


r=^^^^r=vf^ 
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We  then  have 


1      [v 


\/kT( 


TTX)ym 


y/k/(rx)ym-i 


lm 


y/k/{rx)yr 


>/*/(«) 


J/m-1 


where  F(fj)  =  complex  Fresnel  integral  [8] 


TOO  -  ^ 

(X  Jt=0  JT3i 


^m   —  Wm_i  \       /7T 


V^ 


Ay 


m  =  l 


^   I    1/  Vr 

nx 


F  |  ,/ — ym_1 
irx 


(36) 


Consider  now  the  boundary  data  on  the  line  y  =  0.  Assume  that  we  have  a  non- 
uniform grid  0,Xi,x2,-  •  •  ,xp_i/2  —  x.  On  the  interval  (xm_i,a:m),  we  approximate 
the  derivative  as 


9t{t)  = 


um  -  um~l 


Xm         Xm  — \ 

where  um  —  u(xm,0).  At  the  origin  we  have  u°  =  u0.  We  write 


Xp_i/2  gT(T) 


r 

Jr-0 


-At 


\Jxp-\/2  —  r 

We  evaluate  the  first  integral  as 


Jt  =  0       yJXp-i/2  —  T  Jt  =  xt,-i 


I2        9r{r) 


Vxp->/ 


Vxp-i/' 


■dr 


p-3    9r(r)   ^  py  um  -v™-1    n 

JO  Vl    —    T  .    Xm    —    .Tm-I     Jx, 


1 


yX  —  T  m_j  Xm   —  Xm_j   Jxm-i    \JXv_\j2  - 


P-l    ^  _wm-l 


£ 


—  1    •^m  Xm  —  \ 


m  =  l 


Im-l 


=    2E 


P"1    w™  _wm-l 


—  i  Xm         Xrn_j 


m  =  l 


(-2^rp_1/2-r) 

\\Jxv-\/2  —  Xm-\   —  Jlp-l/2  —  Xm  J 


p-l 

=    2E 


um  -  um~1 


=1   \/Xp-l/2  —  Xm-\  +  \/xp-i/2  —  Xr 


(37) 
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Now 


L 


«n-i/a         gT(T) 

Xp-i  y/Xp-i/2  ~ 


.P-1/2  _  „P-1 


zrfr     %     /  — === 

T  Xp-\/2  ~  xp-\   •'a-p-j  yjxp-\j7  ~  T 


Xp- 

1/2  - 

Xp_]   Jxp 

V? 

-1/2  _ 

.UP-1 

(xp 

-1/2  " 

-zP-i)  J 

2  U 

p-1/2 

-  i/P-1 

— 

?2-*p-i   d\ 

7\ 


=    2 

y/x  p-1/2  —  ^p-l 

Substituting  (36),  (37),  and  the  above  in  (35),  we  get 

Zp-l/2 


du 
dy~ 


+ 


v=o 


8j* 


\  7r(a>-i/2  ~  xp-i) 


'p-1/2 


y=0 


Ay 


^2  (u"'  ~  u™-i 


771-1 


F 


■Vm       -  ^ 


T^p-1/2         /  y\J   TT^p-l/2 


■J/m-l 


8;,  g 


um  -  um_1 


7T 


+ 


771  =  1     •v/^Xp-l 


^=7j  +  ^P-l/2-^.-l)  ^ 


8j'A- 


^{xp-l/2  -  Ip-l) 


m""1    (38) 


where 


1  ,  1,  ,A1A 

xp-\/2   —   ^(J'p-1   +  xp)   =3'  xp-l/2  —  -Tp-1    —    n\XP  ~  Xp-l)   —    ^^XP-1 

The  above  equation  is  the  discrete  version  of  a  continuous  boundary  condition  of  the 
form 

—  +  r{x)u  =s{x)  (39) 

oy 


wh 


ere 


r(x)  = 


'8jk  1 


y/X    —    Xp_] 


(40) 


s{x)     = 


n    •     V~^       "771  "771—1 

771=1  •? 


7T2  /  \    V   7TT 


^771     _    ^771-1 


7T      ^   v/l-I,,   +  yfx 


+ 


8j* 


i        \  ft[z  —  xp-\) 


and 


F(x)  =  ["e-JW^d-i 
Jo 


W-1  ,       (41) 


(42) 
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is  the  complex  Fresnel  integral  [8]. 

For  efficient  implementation,  the  PDE  and  the  various  boundary  conditions  will  be 
transformed  into  a  curvilinear  coordinate  system  generated  by  setting  the  lower  ir- 
regular boundary  as  a  constant  curve  curvilinear  coordinate. 
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7.      Transformations  to  a  Curvilinear  Coordinate  System 


Consider  the  narrow  angle  parabolic  equation  with  range  dependent  refractive  index 
(aj  ^  0)  given  in  (18) 

ux  =  — {  a    +  a2—  +  —  \  u         PDE 

(2jA^0-ai)  [  dy       dy7  J 

together  with  the  tropospheric  boundary  condition 

d 


dy 


u(xm,y0)  +  r{xTn)u{xm.y0)  =  s(xm,yQ) 


and  the  impedance  boundary  condition  on  the  irregular  boundary 

uu  +  C\U  =  0 


~7 7 7 7 7 7  7- 


TTT)|:>oS|:> 


^^c      gou^wy 


XoY^kfij^ 


<£/w.CX 


ftrj\*wc/.<*Ay 


We  will  transform  these  to  a  curvilinear  coordinate  system  (£.77) 


c 


^ £ 


4.  c <  / * ,— <: c C-  -yi   -    H 


H >• 


? 


Physical  Domain 


Computational  Domain 
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We  assume  that  we  have  a  transformation  of  the  following  form 

z  =  *(0  ,       y  =  y(Z,v) 

The  various  metrics  needed  in  the  transformed  equations  are  [5] 

011  =  X*  +  J/|  ,  012  =  z*3„  +  y*y„  • 

We  then  have  with  y/g  =  x^yn  -  xvy^  =  x^yn  /  0 

Ui    =    -p{ynui  ~  Viuv)  = — (u(yv-uvys 

1       I  L  \  ^ 

y/g  yn 


l  a 


l 


m  yv^r]  \yj  "  y2v 

Substituting  into  the  PDE  we  get 


Ur,r, Un 

Vr, 


or 


^iVv 


XeCL 


(uiyn  -11^)  = 


1 


uv       y 


r\T, 


nn 


(2j^o-oi)  [  y„       y^  y2v 


zu\ 


■u  + 


02        2/ 


';'/ 


+ 


J/f 


«n  + 


Xi 


{2jk0-a1)"       [\yri        y^  J  (2jk0-a1)    '   yv\  "         {2jk0-a1)y^ 
Letting 

&i     = 


XcQ 


iul 


(2jk0  -  fl! 


62  = 

63  = 


a2 


y 


»7Tj 


1 


+ 


2/£ 


y?  /  (2j*0  -qi)     ^. 


(2j*b-ai)yJ 


we  express  the  PDE  as 


U{  =  biu  +  b2uv  +  63u 


>/r< 


(43) 


The  normal  derivative  on  a  77  =  constant  line  is 

1^(77  =  const.)     =     .  / — uv rig 


v/00 


11 


Ur,    — 


z^r,  +  y&v 


\lx\  +  y? 

^J/t,  -  y^v  yjx\  +  y^x^yrj  -  y^xv) 
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(44) 


Variation  of  an  arbitrary  vector  r  along  the  77  =  const,  coordinate  is 
The  unit  vector  along  the  tangent  on  77  =  const,  is  then 


■;  =  ^£ 


*iX  y^y 


\"l\        y/x\  +  y\       yjx\  +  t/| 


Defin 


nig 


cos0     = 


sin#    = 


\lx\  +  y? 


we  express  (44)  as 


(a:,,  cos#  +  y„  sin  B)u^ 


yTl  cos  0  -  Xr,  sin  0  /x2  +  y2^  cos  5)  _  ^  sin  ^^ 

With  these  substitutions,  the  boundary  condition  at  the  bottom  boundary  uu-\- c^u  = 
0  gets  transformed  to 

(xr,  cos  6  +  yv  sin  0) 

vv ,  u^  +  [yri  cost)  —  xr/  sin  v)cxu  =  U     on  r/  =  0  (4o) 

V*?  +  y\ 

For  jr/  =  0  and  using  J  x2^  +  y%  =  x^/  cosO.  we  get 


u„ sin  #  cos  6u:  +  Ci?/r/  cos  611  =  0     @     r/  =  0 


(46) 


Finally  at  the  top  boundary,  we  have 

uv/yri  +  ru  =  s 


or 


UTl(Zrn'K)  +  riZm-  ■Njyvttm^iZm,  N)   =  yvs((m,N) 


(47) 
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8.      Generation  of  Curvilinear  Coordinate  System 


Consider  a  piece-wise  linear  ground  profile  and  a  horizontal  upper  boundary 


t 


-*■ — < 

— — 1 

r^ 

A' 

t ' 

i — « — 

6" 

1 — , 

, — t — *, 
C' 

1 « 

1 — *- 
0' 

1 

X 

»              < 

1 

/eX. 

^r^C 

)VNs 

/ — t 

k 

c 

Non- Rectangular, 
Uniform  Mesh 

Physical  Domain 


u. — i.  * 


7\ 


-N 


r~* " — r c — f  / 


11 


1 


'r^rr 


A-r 


"Tl 


AYj 


7 7 w—7 7- 

8  c 


fc1  = 


-0 


D 


Rectangular, 
Uniform  Mesh 

Computational  Domain 


»'      f      / 


*t>x, 


xt>,.X 


^-—Z L 1_ 


*l+l  '^1+1 


•-7 t — 7 r 


We  generate  the  (x,  y)  coordinates  of  an  interior  point  by  using  a  linear  interpolation. 
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Letting  Ax,  =  (x,+i  -  x,),  Ay,  =  (yt+]  -  y,),  and  A£,  =  (&+1  -  &),  we  write 

Ax, 


x  =  x,  + 


y  =    i 


A6 
1. 

Ar 


(£  "  6) 


e.-<e<e.-+i 

-M 

0<t]  <N 

(48) 


At  any  interior  point,  £,  <  £  <  f,+1,  the  various  metrics  are  evaluated  as 

Ax, 


xt     = 


A6 


0 


(49) 


2/4 


-* 


v  \  Ay, 

Ae, 


yr,  = 


1 

A7 


j/o-  y. 


Ay, 

A£, 


«-&: 


y™  =  o 


At  the  boundary  points,  we  use  the  central  difference  formulas  to  arrive  at 

xI+2  —  xi       Ax,  +  1  +  Ax, 


*d(  =  Zi) 


£+2-£,  '  '  A(,+]  +  A£ 


t]  \  Ay,+1  +  Ay, 
ydt=Zi)    =       !-T7 


(50) 


(51) 


AV  A£l+1  +  A£, 

Note  that  the  analytical  expressions  yield  a  discontinuous  value  for  these  derivatives 
at  the  boundary  points.  We  use  the  analytical  expressions  only  to  generate  grid 
points  and  use  the  central  difference  formulas  to  arrive  at  the  derivatives  w.r.t.  f. 
In  other  words,  once  the  grid  points  are  generated  on  the  lines  AA' .  BB' ,  .  . . ,  etc., 
we  assume  that  the  space  is  smoothly  connected  through  the  grid  points.  In  the 
numerical  implementation  using  Crank-Xicolson  implicit  scheme  [6].  the  metrics  are 
needed  at  the  midpoint  w.r.t.  £.  i.e.,  at  £  =  £,  +  A£t/2  and  the  interior  point  formulas 
are  applicable.  For  a  uniform  mesh  in  the  computational  domain,  A£,  =  1,  r\  =  9, 
9  =  0,1, 2,...,iV 


x( 

— 

Ax, 

Vi 

= 

i  / 

yr, 

= 

aH90" 

y 

= 

('-« 

= 

Vi+i  +  y. 

2 

- 

y(?) 

=  yrr 

Ay, 


yi+i  +  y* 


y^  +  i)-y*(9)  =  - 


Ay, 
N 


(52) 
(53) 

(54) 


y.+i  +  y> 


+  £* 


+  gyr,  =  yo  +  (9  -  Ar)yr, 


(55) 
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9.     Numerical  Implementation  by  Crank-Nicolson  Scheme 

Consider  the  narrow  angle  PE  together  with  the  boundary  condition 

u<  =  bxu  +  62u,  +  b3um         (PDE) 

uri  —  {  —  sin  0cos  0  }  U(  +  Ciy^  cosOu  =  0     at  n  =  0         (BC2) 

»,  +  rj/,11  =  yri  s     at  n  =  A'         (BCi ) 
We  would  like  to  implement  the  above  using  a  Crank-Nicolson  implicit  scheme  [6] 


Av\-  I 


p-l 


A?  = 


f>-i 


f> 


t 


% 


V 


We  use  the  notation  up  =  u(^p.T]q)  —  u(xp.yq) 
The  various  derivatives  assuming  A£  =  Ar/  =  1  are 

U^p-1/2.T]q)       =      UPq-VPq'1 

Un(p--.q)     =     u^p-i/2,1,)  ~  2  M^p-^Vg)  +  uv(ZP'Vc,)} 


or 


1  r 


,p-i       „,p-i 


U*q+1  -Uj_!   +W 


9+1  U9-l 


W??  V       2'  V  =  4  ^+]1  +  U'+1  ~  (U'-jl  +  U'_1) 

p  -  i?)  =  £  k+i  -  2< + <-i  +  «*s  -  2<l + <=! 
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Substituting  into  the  PDE,  we  get 

1     \  t£+  up~l 


v  v— 1 

ir  —  ir 


61  (^  ~  2"9 


9  V       2    V  2 

1 


^2  [p-  -.?)   , 


Rearranging  the  terms  we  get 


K^*.K«-  ,+fc-y  <+5  fc-yK« 


1  /  62  •  !  ^ 

2  V  2 


+i(?H|tC!+  i-^  +  tK'1  +  ;I*»-?)«P!!S0 


Now  let 


1  /       1    v_1/2 

»-M). 

We  then  have  for  p  =  1,2,  •  ■  •,  q  =  1,2,  •  •  ■,  N  -  1 

auj+1  -  (1  +  0)u;  +  7<-i  =  "«<;]  -  (1  -  /?K_1  "  7<:J  (56) 

However,  we  will  extend  the  applicability  of  this  equation  over  q  —  0,1,-  ■  •  Ar  to 
accomodate  the  derivative  boundary  conditions  on  the  lower  and  upper  boundaries. 
For  q  —  0.  we  have 

au\  -  (1  +  /?K  +  7uP-i  =  -our1  "  (1  "  ^uS"1  -  7"P-~11  (57) 

for  q  =  N,  we  have 

a<v+1  -  (1  +  $)u'N  +  7^-1  -  -<i  -  (1  -  Wiijr1  -  7<"-1i  (58) 
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From  the  BC2  at  the  lower  boundary  we  have 

1 


p-1/2 


sin  $  cos  0  1  (up0-up0~l) 


.    ,  mP-1/2  (V^±U^1\         n 

+  {ciVv  cos  0)o         ( 2 J 


Multiplying  this  with  4(7)0"       and  adding  to  (57) 


(a  +  7K-     1  +  /?  -  27 


cj  yn  cos  0  —  2—  sin  0  cos  6 
xi 


up0 


Letting 


=  -(q  +  7K      "     1-/5  +  27 


O        =      Q  +  -) 


c^  cos  0  +  2—  sin  0  cos  6 
xt 


,p-i 


/?'    -    /?-2Oy7?cos0(c1- 
/?"    =    /?-2<yy„cos0(Cl  + 


2sin6>' 
2sinfl' 


we  ma}'  write  the  above  equation  as 

a'u\  -  (1  +  /?X  =  -o'til"1  "  (1  "  0'X"1 
From  the  BC^  at  the  top  boundary  we  have 


'  <»<    <    k 


-*—£- 


P-t 


,    r-    N 
P 


(59) 
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2  12  2f" 


16  jk    uA-  +  i/A' 


p-1 


\J  7rAxp_i 


^V    m  =  l 


F 


.N 


Tr^p-1/2 


•2/m       ~  F 


^  tt^p-i/: 


2/m-l 


1/  XJ  -  U 


771-1 


N       u  A" 


+ 


16;* 


'8J*  '" '  


p-1 


or 


j& 


i  /  *> 

multiplying  with  —  aN        and  adding  to  (58) 

i&  + (7 +  «)«&-, 


/2 


1  +  0  +  y„8a 
-  /?  -  y„8a 


jk 


jk 


\  TTAXp.! 

uAT1-(7  +  a)UAr_11-4ay^-1/2 


Letting 


\J  7rAxp_! 


A    =     .3  +  2^80. 


J* 


^A 


Xp-] 


1      =     7  +Q 
we  rewrite  the  above  equation  as 

(1  +  xyN  -  y'u^  =  (1  -  AKr1  +  y^ri  +  40^-1/2 

where 

Ay  , 


(60) 


5 


p-1/2 


X]    (Um   ~  Wm-l) 


F 


A 


* 


7TXp_i/2 


■J/m  I    -  F 


\  KXp-i/: 


■Vm-l 


m         „.rn~  1 


UAr    -    «A, 


'       ^      m-\    {y/Xp-1/2  —  ^m   +  \/xp-\/2  —  ^m-lj 


+ 


\ 


7rAxp_! 


p-i 
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an( 


Ax 


v-\ 


xp-\/2 


■i  +  2Alp~l 


'\Xp  "T  3-p— 1  ) 


We  augment  the  equations  given  in  (56)  for  q  =  1, 2, . . . ,  JV  —  1  with  (59)  and  (60)  for 
q  =  0  and  A',  respectively  to  define  for  q  =  0. 1,  2, . . . ,  N .  The  system  of  equations  so 
defined  can  be  expressed  as  a  matrix  equation  of  the  form 


X    X 
X    X    X 
X    X    X 


X    X 


X    X 
X    X    X 
XXX 


X    X 


-1  1 


,p-l 


,p-l 

lN 


+ 


4a2/„s 


p-l/2 


(61) 


where  X  denotes  a  non-zero  entry.  The  tridiagonal  matrix  on  the  left  hand  side  of 
(61)  can  be  inverted  efficiently  to  yield  a  solution  on  line  £  =  £p  in  terms  of  the  field 
values  on  the  line  £  =  £p_i.  Equation  (61)  can  be  used  to  march  forward  in  range 
starting  from  initial  data  specified  on  £  =  0. 


34 


References 


[1.]  F.  B.  Jensen,  et  al.,  Computational  Ocean  Acoustics,  Chapter  6,  New  York:  AIP 
Press,  1993. 

[2.]  E.  Zauderer,  Partial  Differential  Equations  of  Applied  Mathematics,  New  York: 
John  Wiley  k  Sons,  1983. 

[3.]  T.  B.  A.  Senior,  "Impedance  boundary  conditions  for  imperfectly  conducting 
surfaces,"  Appl.  Sa.  Res.,  Sec.  B.  Vol.  8,  pp.  481-436,  1960. 

[4.]  S.  W.  Marcus,  "A  hybrid  (finite  difference  surface  Green's  function)  method  for 
computing  transmission  losses  in  an  inhomogeneous  atmosphere  over  irregular 
terrain,"  IEEE  Trans.  Antennas  Propagat..  Vol.  40,  No.  12.  pp.  1451-1458. 
Dec.  1992. 

[5.]  J.  F.  Thompson,  Z.  U.  A.  Warsi,  and  C.  W.  Mastin,  Numerical  Grid  Generation. 
New  York:  North  Holland,  1985. 

[6.]  G.  D.  Smith,  Numerical  Solution  of  Partial  Differential  Equations:  Finite  Dif- 
ference Methods,  Oxford:  Clarendon  Press,  1985. 

[7.]  C.  A.  Balanis.  Advanced  Engineering  Electromagnetics,  New  York:  John  Wiley 
k  Sons,  1989. 

[8.]  M.  Abramowitz  and  I.  A.  Stegun,  Handbook  of  Mathematical  Functions,  New 
York:  Dover  Publications,  1972. 


35 


INITIAL  DISTRIBUTION  LIST 

No.  Copies 

1.  Defense  Technical  Information  Center  2 
Cameron  Station 

Alexandria,  VA  22314-6145 

2.  Dudley  Knox  Library,  Code  52  2 
Naval  Postgraduate  School 

Monterey,  CA  93943-5002 

3.  Chairman,  Code  EC  1 
Department  of  Electrical  and  Computer  Engineering 

Naval  Postgraduate  School 
833  Dyer  Road,  Room  437 
Monterey,  CA  93943-5121 

4.  Professor  Ramakrishna  Janaswamy,  Code  EC/Js  8 
Department  of  Electrical  and  Computer  Engineering 

Naval  Postgraduate  School 
Monterey,  CA  93943-5121 


36 


DUDLEY  KNOX  LIBRARY 


3  2768  00331658  9 


